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HIGHLIGHTS 


•  A  mathematical  model  for  the  thermal  behavior  of  a  wound  cylindrical  Li-ion  cell  was  developed. 

•  The  new  numerical  method  significantly  saves  computation  time  and  memory  allocation. 

•  The  coordinate  transform  and  variable  extrusion  algorithms  were  validated. 
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A  new  computational  method  is  proposed  that  can  be  used  to  reduce  the  numerical  difficulties  in 
modeling  the  electrical  and  thermal  behavior  of  a  spirally  wound  Li-ion  cell.  By  analyzing  the  winding 
locus  of  the  electrodes,  some  important  geometric  relationships  of  the  spiral  surfaces  are  identified,  and 
algorithms  for  coordinate  transform  and  variable  extrusion  between  2-D  and  3-D  domains  are  derived. 
Our  method  reduces  the  computation  time  and  memory  requirements  needed  to  simulate  the  cell 
performance.  The  accuracy  of  our  method  was  validated  by  model-to-model  comparisons. 

©  2013  Published  by  Elsevier  B.V. 


1.  Introduction 

Today,  most  of  the  safety  issues  for  Li-ion  cells  and  batteries  are 
related  to  the  electrical  and  thermal  properties  of  the  composite 
electrodes.  [1—4]  Therefore,  mathematical  models  that  deal  with 
the  charge  balance  and  heat  transfer  in  cells  have  become 
important  tools  for  cell  design.  However,  the  implementation  of 
these  models  is  significantly  affected  by  the  cell  geometry  and 
electrode  configurations.  There  are  basically  two  types  of  electrode 
structures  in  Li-ion  cells,  the  laminated  or  pouch  structure  and  the 
wound  structure.  For  the  laminated  electrode  structure,  all  the 
electrodes  are  flat  and  arrayed  in  parallel,  and  the  transfer  of 
charge  and  energy  goes  along  the  orthogonal  directions  in  the 
Cartesian  coordinate  system.  Therefore,  building  and  meshing  the 
computational  domains  for  laminated  electrodes  can  be  done 
easily  in  commercial  CAD  software.  For  the  wound  electrode 
structure,  however,  the  electrodes  are  put  through  a  complicated 
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winding  locus,  and  the  electrical  potentials  of  the  solid  phases  are 
non-continuous  along  the  normal  direction  of  the  winding  surface 
due  to  the  insulating  separator  layer.  Therefore,  electrical  insu¬ 
lation  zones  with  full  winding  details  must  be  included  in  the 
model  domains  for  wound  electrodes,  and  the  meshing  of  these 
domains  becomes  extremely  difficult  if  the  electrodes  have  many 
winding  turns. 

Due  to  the  numerical  difficulties  discussed  above,  most  3-D 
physical  models  have  been  developed  for  cells  with  laminated 
electrodes  [5-9],  and  modeling  for  the  wound-structure  elec¬ 
trodes  have  been  rarely  reported.  Gomadam  et  al.  [10  developed 
a  2-D  thermal  model  for  a  cylindrical  jell-roll  based  on  the  spiral 
geometry  relations,  but  the  charge  balances  in  the  current  col¬ 
lector  foils,  which  are  supposed  to  be  the  most  difficult  parts  for 
the  wound  electrode  models,  were  neglected.  Santhanagopalan 
et  al.  [11]  and  Ye  et  al.  [12]  also  developed  physics-based  models 
for  the  wound  cells,  but  their  models  are  limited  to  2-D.  In  a 
recent  paper  by  Lee  et  al.  [13],  a  multi-scale  multi-dimensional 
(MSMD)  model  for  wound  cells  was  reported;  however,  their 
approach  involves  two  mesh  groups  in  which  the  size  for  ge¬ 
ometry  units  are  no  larger  than  twice  the  electrode  pair 
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thickness,  and  thus  requires  a  large  number  of  mesh  elements 
which  would  make  the  model  computationally  expensive. 

In  this  work,  a  new  modeling  approach  different  from  all  the 
above-mentioned  ones  is  developed.  The  charge  balance  and 
heat  transfer  in  a  jelly-roll  are  solved  separately  in  different 
computational  domains,  and  coupled  through  a  coordinate 
transform  and  variable  extrusion  algorithm  derived  from 
analytical  geometry  relationships.  The  mesh  requirement  for  the 
3-D  domain  is  significantly  lowered  by  using  these  geometric 
algorithms. 

2.  Mathematical  approach 

2.1.  Basic  geometric  relations 

Our  model  was  developed  from  an  assumption  that  all  the  cell 
sandwich  layers  (current  collectors,  coatings,  and  the  separator)  in 
the  jelly- roll  of  al 8650-type  cylindrical  cell  are  of  uniform  thick¬ 
nesses  and  are  wound  through  an  Archimedean  spiral  locus  which 
can  be  described  by  the  following  polar  coordinate  equation: 

r-r.+^O  (,) 

where  r  is  the  radius  coordinate,  6  is  the  angular  coordinate,  ro  is 
the  initial  radius  of  the  spiral  at  0  =  0,  and  D  is  the  separation 
distance.  The  geometry  details  for  a  spiral  surface  are  presented 
in  Fig.  1,  where  the  surface  is  wound  spirally  around  the  z  axis 
and  the  cross  section  on  the  x-y  plane  is  an  Archimedean  spiral 
curve  described  by  Equation  (1).  The  Archimedean  spiral  has  the 
property  that  any  ray  from  the  origin  intersects  successive 
turnings  of  the  spiral  in  points  with  a  constant  separation  dis¬ 
tance  D,  and  the  difference  in  the  polar  angle  between  two 
intersecting  points,  A0,  is  2tc.  In  Fig.  1,  f  is  the  tangential  arc 
length  along  the  spiral  curve  and  differentiation  of  ?  can  be 
expressed  as  follow: 


A  functional  relationship  between  the  arc  length  ?  and  the  polar 
angle  0  can  be  derived  by  integrating  Equation  (2): 


(3) 

Equation  (3)  sows  that  f  is  a  unique  function  of  0;  therefore,  at  a 
fixed  z,  any  variable  defined  on  the  spiral  surface  can  be  expressed 
as  a  function  of  the  polar  angle  0.  Due  to  the  small  thicknesses  of  the 
Li-ion  cell  electrodes,  it  can  be  found  that  (r0  +  (D/2tu)0)2  is  at  least 
three  orders  larger  than  (D/2tu)2  for  all  0,  and  Equation  (3)  can  be 
approximated  as: 


therefore,  an  approximate  functional  relationship  between  0  and  ? 
can  be  derived  from  Equation  (4): 


2tz 

As  shown  in  Equation  (5),  at  a  fixed  z,  any  variable  defined  in  the 
spiral  surface  can  also  be  expressed  as  a  function  of  the  arc  length  f. 

According  to  references  [5-9],  the  current  collector  foils  can  be 
regarded  as  surfaces  due  to  the  small  thickness,  and  the  electrical 
potential  of  a  current  collector  distributes  only  in  the  two  tangen¬ 
tial  directions  of  the  surface  (?  and  z  directions).  Therefore,  the 
charge  balance  equations  on  the  current  collector  surfaces  are: 


where  4>j  is  the  electric  potential  of  the  current  collector  of  elec¬ 
trode  j,  Gj  is  the  electric  conductivity  for  the  current  collector  of 
electrode  j,  and  ivj  is  the  volumetric  current  source  in  the  current 
collector  of  electrode  j.  It  can  be  found  that  4>j  solved  from  Equation 
(6)  is  a  function  of  ?  and  z;  therefore  at  a  fixed  z,  $j  can  be  expressed 
as  a  function  of  the  polar  angle  0. 

As  discussed  later  in  our  work,  certain  electrical  and  thermal 
variables  need  to  be  evaluated  along  the  spiral  surface  and 
extruded  into  the  3-D  Cartesian  coordinate  system.  The  extru¬ 
sion  of  variables  can  be  performed  through  a  coordinate  trans¬ 
form  algorithm.  As  shown  in  Fig.  1,  t|/(0)  is  a  variable  defined  on 
the  spiral  surface  at  a  fixed  z,  and  (x,y)  is  a  specific  Cartesian 
coordinate  point  which  has  the  same  z  coordinate  as  <j/(0).  Point 
(x,y)  could  be  located  either  on  or  outside  the  spiral  surface.  The 
first  step  for  variable  extrusion  is  to  project  point  (x,y)  to  the 
spiral  surface  along  a  specific  direction  defined  by  a  ray  from  the 
origin  and  passing  (x,y),  the  elevation  angle  0e  (where 
0  <0e  <  2tt)  for  the  projecting  direction  can  be  calculated  as 
follow: 


Fig.  1.  Schematic  for  a  surface  wound  through  an  Archimedean  spiral  locus. 


222 


M.  Guo,  R.E.  White  /  Journal  of  Power  Sources  250  (2014)  220-235 


arccos  (  *  ,)  if  y  >  0 

2iz  -  arccos  (  .  =  „  )  if  v  <  0 

W*2+y2/ 

The  direction  of  projection  is  illustrated  by  the  red  (in  the  web 
version)  arrow  in  Fig.  1  which  can  be  approximately  regarded  as 
normal  to  the  spiral  surface.  Next,  define  an  integer  nw  which  de¬ 
fines  an  index  for  the  nearest  turning  of  the  spiral  surface  to  point 

(xoO: 


nw(x,y) 


v^TF-ro-^E 

D 


(8) 


where  the  square  brackets  “[]”  in  this  equation  denote  an  operator 
to  find  the  nearest  integer  to  a  real  value.  Thus,  the  polar  angle  for 
the  projection  of  point  (x,y)  on  the  spiral  surface  can  be  calculated 
as: 


Coating  of  positive  electrode 
Coating  of  negative  electrode 
Separator 

-  Current  collector  of  positive  electrode 

-  Current  collector  of  negative  electrode 


<t>.(0-2x)  « 


®  (»)  ■ 


<M«) 


<t>_(0  +  2n) 


Stack  (9-2 n,9)  Stack  (9,9)  Stack  (9,9  +  27t) 


6(x,y)  =  2  7t-nw(x,y)  +  BE(x,y)  (9) 

Using  Equation  (9),  a  transform  between  the  polar  and  Cartesian 
coordinates  can  be  established: 

W)  =  t'my))  =  ^(x,y)  (10) 

The  algorithm  described  by  Equations  (7)— (10)  shows  that  any 
variable  evaluated  on  a  spiral  surface  can  be  extruded  to  a  Cartesian 
coordinate  position  in  the  space  enclosed  by  the  spiral  surface. 

The  charge  balance  Equation  (6)  is  defined  on  a  3-D  spiral  sur¬ 
face,  but  the  dependent  variables  can  be  continuously  distributed 
in  the  two  tangential  directions  of  the  surface.  Therefore,  the 
computational  domain  for  Equation  (6)  can  be  transformed  into  a 
2-D  flat  surface  (see  Fig.  2),  the  f-z  plane;  the  solutions  for  the 
charge  balance  are  approximately  same  on  both  the  spiral  and  the 
planar  surfaces.  This  assertion  will  be  validated  later  in  this  article. 

2.2.  The  numerical  difficulties  for  conventional  modeling 
approaches 

The  configurations  for  a  sector  of  the  spirally-wound  electrode 
pair  in  a  jelly-roll  are  presented  in  Fig.  3,  where  the  region  inside 
the  dashed  box  shows  the  array  of  layers  intersecting  a  ray  from  the 
origin.  As  shown  in  Fig.  3,  the  parts  between  two  adjacent  opposite 
current  collectors  are  defined  as  cell  stacks  and  each  cell  stack  in¬ 
cludes  a  positive  electrode  coating  layer,  a  separator  layer,  and  a 
negative  electrode  coating  layer.  The  electrical  potentials  of  the 
current  collectors  are  expressed  as  functions  of  polar  angles;  and 
the  phrase  “Stack(0,0')”  stands  for  a  cell  stack  in  which  the  polar 
angles  for  the  positive  and  negative  electrode  current  collectors  are 
respectively  6  and  O’  (where  O’  =  0,  6  +  27tj,  or  6  -  2iz).  The  symbol  “ 
i"n  (0  o')  "  denotes  the  transversal  current  density  through  the  coating 
layers  of  stack  (0,0'). 


Fig.  3.  Configuration  for  wound  electrode  pair. 


Fig.  4.  Mesh  pattern  for  the  computation  domains  of  jelly-roll  where  electrical  insu¬ 
lation  zones  are  included. 


Spiral  surface  Planar  surface 

Fig.  2.  Transform  of  computational  domains  for  electrical  equation  for  electrode  j. 


Fig.  3  shows  that  between  two  neighboring  stacks,  one  of  the 
two  current  collector  potentials,  $+  or  must  be  discontinuous; 
for  example,  from  Stack(0  -  2tu,0)  to  Stack(0,0),  the  polar  angle  of  &+ 
jumps  by  2tt  and  therefore  T>+  is  discontinuous  through  these  two 
stacks.  However,  the  commonly-used  numerical  methods  (finite 
element  method,  finite  volume  method,  etc.)  require  all  variables  to 
be  continuous  in  a  mesh  element.  Therefore,  if  the  charge  balance 
Equation  (6)  is  directly  solved  on  the  3-D  spiral  surface,  the 
maximum  radial  dimension  for  each  mesh  element  must  not 
exceed  the  thickness  of  a  cell  stack.  Generally,  the  thickness  of  a  cell 
stack  is  smaller  than  the  dimensions  of  a  jelly-roll  by  two  orders  of 
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Tab-contact  surface 


(a)  Configuration  of  individual  electrodes 


Negative  electrode 


Fig.  5.  The  configurations  for  (a)  individual  electrodes,  (b)  electrode  pair. 


magnitude,  and  a  very  large  number  of  mesh  elements  are  required 
to  ensure  sufficient  accuracy.  Fig.  4  shows  the  mesh  pattern  for  the 
modeling  domain  of  al8650-type  jelly-roll  generated  by  COMSOL 
4.3b.  In  the  jelly-roll  domain  presented  in  Fig.  4,  two  20-turn  spiral 
surfaces  are  built  as  the  electrical  insulation  interior  boundaries, 
and  a  total  of  511,241  domain  elements  are  included  to  meet  the 
default  lowest  mesh  requirements  by  COMSOL,  and  this  extremely 
fine  mesh  pattern  would  make  the  model  computationally 
expensive. 

2.3.  Solutions  for  computation  problems 

Based  on  the  ideas  discussed  in  Section  2.1,  we  propose  a  new 
method  which  can  greatly  reduce  the  computational  load  for  a  3-D 
spirally  wound  cell  model.  The  electrical  equation  (6)  is  solved  on  a 
transformed  2-D  planar  domain  for  an  electrode  pair  which  is 
expressed  as  the  f-z  plane.  The  Joule  heat  source  and  the  electro¬ 
chemical  heat  source  are  also  evaluated  on  the  2-D  electrode  pair 
domain  and  expressed  as  functions  of  6  and  z;  and  these  heat  source 
terms  are  extruded  to  the  jelly-roll  domain  built  in  the  3-D  Cartesian 
coordinate  system  to  be  coupled  with  the  3-D  thermal  model. 

2.4.  The  electrode  configurations 

Configurations  for  the  two  individual  electrodes  are  presented 
in  Fig.  5(a).  Each  electrode  includes  three  different  coating  regions, 


the  bare  foil  region  where  no  coating  is  attached  to  the  current 
collector  foil,  the  single-sided  region  where  only  one  side  of  the 
current  collector  has  an  active  material  coating,  and  the  double¬ 
sided  region  where  both  sides  of  the  current  collector  are  coated. 
In  the  bare-foil  region,  there  is  a  zone  called  the  tab-contact  surface 
on  which  the  current  enters/leaves  the  current  collector  form/to  the 
external  tabs.  As  shown  in  Fig.  5(b),  by  overlapping  the  coating 
regions  of  the  two  electrodes,  the  2-D  domain  for  electrode  pair  is 
obtained.  In  Fig.  5(b),  the  NS-PD  region  is  the  matching  region  for 
the  single-sided  negative  electrode  and  the  double-sided  positive 
electrode,  the  ND-PD  region  is  the  matching  region  for  the  double¬ 
sided  negative  electrode  and  the  double-sided  positive  electrode, 
and  the  ND-PS  region  is  the  matching  region  for  the  double-sided 
negative  electrode  and  the  single-sided  positive  electrode.  The 
top-views  for  the  matching  of  the  coating  regions  in  an  electrode 
pair  are  presented  in  Fig.  6. 

2.5.  Electrochemical  model  for  cell  stacks 

In  this  paper  we  have  chosen  to  use  an  empirical  equivalent 
circuit  model  [14]  (ECM)  to  simulate  the  potential  and  current  re¬ 
sponses  of  each  cell  stack.  According  to  reference  14,  the  equivalent 
circuit  (see  Fig.  7)  includes  two  R—C  branches  {Rp  -  Cp  and  Rn  -  Cn), 
a  series  resistor  (Rs),  and  an  open  circuit  voltage  source  (U).  The 
governing  equations  and  validation  of  the  ECM  are  discussed  in 
Appendix  A. 
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(a)  Top- view  for  individual  electrodes 


End  of  the  single-sided  region  and 


for  positive  electrode 


(b)  The  matching  of  coating  layers  in  the  electrode  pair 


Fig.  6.  Top-views  for  the  matching  of  coating  regions  in  an  electrode  pair. 


2.6.  The  current  source  terms 

According  to  references  [5—9],  the  volumetric  current  source  in 
the  current  collectors, iv,j  (j  =  +,-),  are  expressed  by  the  transversal 
current  density  scaled  by  the  thickness  of  current  collector. 
Therefore,  iVj  has  a  different  expression  in  the  different  coating 
regions  of  the  2-D  electrode  pair  domain  presented  in  Fig.  5(b).  The 
arrays  of  electrode  component  layers  for  the  NS-PD,  ND-PD,  and 
ND-PS  regions  of  a  wound  electrode  pair  are  respectively  presented 
in  Fig.  8(a)  through  (c).  Note  that  to  make  the  illustrative  sche¬ 
matics  clear;  the  wound  electrode  pairs  in  Fig.  8  are  presented  as 
one-turn  spirals,  however,  the  practical  jell-roll  often  has  many 
more  winding  turns.  In  the  NS-PD  region  shown  in  Fig.  8(a),  the 
negative  electrode  current  collector  has  only  one  side  with  the 
outward  transversal  current  density  and  the  positive  electrode 
current  collector  has  both  sides  with  the  inward  transversal  current 


Stack  [0,6') 

l 

o_(0')  O+(0) 


lN,(0,0') 


Cp  Cn 


Coating  of  positive  electrode 
Coating  of  negative  electrode 
Separator 

Current  collector  of  positive  electrode 
Current  collector  of  negative  electrode 


Fig.  7.  Equivalent  circuit  for  a  cell  stack. 
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(a)  Electrode  component  arrays  for  the  NS-PD  region 


(b)  Electrode  component  arrays  for  the  ND-PD  region 


Coating  of  positive  electrode 
Coating  of  negative  electrode 
Separator 

Current  collector  of  positive  electrode 
Current  collector  of  negative  electrode 


- 1 - 1 - 

Stack  ( 0,0+2it )  Stack  (0 +2it  ,0 +2n) 

0+  (o)  O_(0  +  2;r)  0+  (0+2n ) 


lU,(0,O+2x) 


lK,(O+2z,0+2x) 


""  /?  **P  Kn  ft 


(c)  Electrode  component  arrays  for  the  ND-PS  region 


Fig.  8.  The  array  patterns  for  electrode  component  layers  in  different  coating  regions:  (a)  for  NS-PD  region,  (b)  for  ND-PD  region,  (c)  for  ND-PS  region. 


density;  therefore,  the  volumetric  current  source  terms  at  the  polar 
angle  0\  =  6  -  2tt  are  expressed  as: 


iv,+(0  i) 


IN, (0i  ,00  +  IN,(01,01+2tu) 


(11) 


shown  in  Fig.  8(b),  the  negative  electrode  current  collector  has  both 
sides  with  the  outward  transversal  current  density  and  the  positive 
electrode  current  collector  has  both  sides  with  the  inward  trans¬ 
versal  current  density;  therefore,  the  volumetric  current  source 
terms  at  the  polar  angle  6  are  expressed  as: 


iv,-(*  i) 


’n.qMi) 


<5_ 


(12) 


iv,+(0) 


lN,(0,0)  +  lN,(0,0+27C) 

K 


(13) 


where  d+  and  o_  are  respectively  the  thicknesses  for  the  current 
collectors  of  positive  and  negative  electrodes.  In  the  ND-PD  region 


>v,-(0)  = 


lN, (0-271,0)  +  1'n,(0,0) 

T- 


(14) 
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Rp  Kn 
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Fig.  9.  Definition  of  the  average  heat  generation  rate  for  two  neighboring  cell  stacks. 


In  the  ND-PS  region  shown  in  Fig.  8(c),  the  negative  electrode 
current  collector  has  both  sides  with  the  outward  transversal  cur¬ 
rent  density  and  the  positive  electrode  current  collector  has  only 
one  side  with  the  inward  transversal  current  density;  therefore,  the 
volumetric  current  source  terms  at  the  polar  angle  82  =  8  +  2n  are 
expressed  as: 


iv,+  (02) 


*n,(M2) 

<5+ 


(15) 


iv,-(02) 


*'n,(#2-2tu,02)  +  lN,(02,02) 

8l 


(16) 


In  the  tab-contact  surfaces  of  the  bare-foil  regions,  the  current 
source  terms  are  expressed  as: 


iv,+ 


*App 

^+<5+ 


(17) 


2.7.  The  heat  source  terms 


According  to  the  equivalent  circuit  model,  the  volume-average 
electrochemical  heat  through  the  thickness  of  the  cell  stack  in 
Fig.  7  can  be  calculated  as: 


QeC ,(0/)  ~ 


^)  [*+(*)- MO  ~»] 


(19) 


where  5S  is  the  thickness  for  a  cell  stack.  The  volume-average  Joule 
heat  through  the  cell  stack  is  expressed  as: 


0.5<5+<7+ 


~  —5, 


fd<P+ 
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/0<2>+ 
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\  9? 

J  +' 

(  0 Z 
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0.5 

5S 
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(20) 


and  the  total  heat  generation  rate  throughout  cell  stack  (6,6')  is: 


iv’~  =  (18) 

where  A+  and  A_  are  the  area  of  the  tab-contact  surfaces  for  the 
positive  and  negative  electrodes,  and  Japp  is  the  applied  current  on 
the  cell  which  is  defined  as  positive  for  charge  and  negative  for 
discharge.  The  current  source  is  zero  at  other  positions  in  the  bare- 
foil  regions. 


Qh,(0/)  -  Qec,(0/)+Qj,(0/)  (21) 

The  average  heat  generation  rate  (see  Fig.  9)  for  the  two  adjacent 
cell  stacks  which  share  the  same  positive  electrode  current  col¬ 
lector  is  obtained  as  a  unique  function  of  polar  angle  8: 

Qh  (6)  =  Qh.(^)  +2Qh  W)  (22) 


Table  1 

Expressions  for  volumetric  current  source  and  heat  generation  rates  in  different  coating  regions. 


Coating  regions 

Current  source  terms 

Heat  source  terms 

Bare-foil  region  for  negative  electrode 

=  (tab-contact  -  surface) 

Qh(0)  =  v-  [(irU  + 

(^L)21 

iy.-(O) 

=  0  (other  locations) 

_  Oh,(0,9) +Qh,(9,0+27t) 

NS-PD  region 

iv,+(«) 

+iN,(V+2TC) 

(5+ 

“  d~ 

_  Qh,(0,0)  +Qh,(0,0+27r) 

ND-PD  region 

iv,+(0) 

1’n,(0,0)  +iN,(0,0+2ir) 

“  <5+ 

iv,-(0) 

lN,(6- 2%,0)+lN,(8,0) 

6 

ND-PS  region 

iv.+(«) 

_  *N,(9.«) 

Qh(0)  =  Qh,(d,d) 

IN,(0-2tc,0)+IN,(0,0) 

6 

Bare-foil  region  for  positive  electrode 

iv,+(«) 

=  (tab-contact  -  surface) 

QH(0)-ff+[(t|#)2  + 

0H)2] 

iv,+(«) 

=  0  (other  locations) 
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Note  that  in  the  ND-PS  region,  the  positive  electrode  current 
collector  is  single-sided  and  the  expression  for  the  heat  generation 
rate  is  as  follow: 


QhW  =  Qh  m  (23) 

In  the  bare-foil  regions,  the  electrochemical  heat  is  zero  and  the 
heat  generation  rates  only  include  the  Joule  heat;  therefore  for  the 
positive  electrode 


Qh(0) 


and  for  the  negative  electrode 


Qh  (0) 


21 


(24) 


(25) 


The  expressions  for  the  volumetric  current  sources  and  heat 
generation  rates  in  different  coating  regions  are  summarized  in 

Table  1. 


2.8.  Thermal  model 

Since  all  of  the  electrode  component  layers  are  thermally 
conductive,  the  heat  fluxes  in  the  jelly- roll  are  continuous  in  all 
three  directions;  therefore,  the  thermal  model  must  be  solved  in  a 
3-D  computational  domain.  The  governing  equation  for  the  thermal 
model  is  as  follow: 


+  Qh 


(26) 


where  p  is  the  density  of  the  jelly-roll,  Cp  is  the  specific  heat  ca¬ 
pacity  of  the  jelly-roll,  K  is  the  matrix  for  the  anisotropic  thermal 
conductivity  of  the  jelly- roll,  and  is  the  extruded  volumetric 
heat  source  in  the  3-D  Cartesian  coordinate  system  in  the  jelly- 
roll. 

In  a  jelly- roll,  the  electrode  component  layers  are  arrayed  in 
series  along  the  normal  direction  of  the  spiral  surface  and  in  par¬ 
allel  along  the  tangential  direction  of  the  spiral  surface.  Due  to  such 
an  array  pattern,  the  thermal  conductivity  of  the  jelly-roll  has 
different  values  in  the  tangential  and  normal  directions  of  the  spiral 
surface.  As  shown  in  Fig.  10,  at  a  location  with  polar  angle  6 ,  N 5 
denotes  the  heat  flux  along  the  tangential  direction  of  the  spiral 
surface,  Nr  denotes  the  heat  flux  along  the  normal  direction  from 
the  spiral  surface,  and  the  value  for  6  can  be  determined  from  the 
Cartesian  coordinates  using  Equations  (7)-(9).  The  angle  between 
the  normal  direction  of  the  spiral  surface  and  the  x  axis,  a,  can  be 
derived  as: 

a  =  6  +  arctan(27r^  - ^  (27) 

The  expressions  for  N 5  and  Nr  in  the  Cartesian  coordinate  system 
are  as  follow: 

AT  ,  dT  (  .  dT  dT\ 

Nr  =  =  -hr  l  -  since  —  +  cosa—  (28) 

k  dx  dyj  v  ' 


Fig.  11.  Flowchart  for  the  solution  procedure  using  our  proposed  algorithm. 
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r - l 


(a)  The  3-D  jelly-roll  domain 


M  i  ,  /  or  .  an  f  . 

Nr  =  -fcr-=-kr^°s«-  +  sin«-j  (29) 

where  k ^  and  kr  are  respectively  the  thermal  conductivities  in  the 
tangential  and  normal  directions  of  the  spiral  surface.  The  heat 
fluxes  in  the  x,  y,  and  z  directions  can  be  expressed  as: 


sincdVr  +  cos  aN^ 

~(kT  -  k^j sinacosa^  -  (/<rsin2a:  +  k^c os2 a'j  ^ 


(31) 


Nx  =  coscdVr  -  sin aN^ 


-  (l<vcos2a  +  k^sin2^  ^  -  (kT  -  k^j 


- /<p  )  sinacosa  ^  (30) 


Nz  =  -k. 


a  t 

dz 


(32) 


Equations  (30)— (32)  can  be  expressed  in  vector-matrix  format: 
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Fig.  13.  Mesh  pattern  for  the  3-D  jelly-roli  domain. 


Table  2 

Values  of  parameters.3 


Symbol 

Value 

Unit 

C 

1.17 

Ah 

CP 

750 

J  K1  kg-1 

D 

3.25  x  10"4 

m 

ic 

15 

Am"2 

kT 

0.8 

W  K1  nr1 

h 

30 

W  K1  nr1 

Nw 

20 

r0 

2  x  10"3 

m 

<5+ 

20  x  10"6 

m 

<5_ 

15  x  10~6 

m 

1.625  x  10"4 

m 

P 

2720 

kg  ITT3 

a+ 

3.8  x  107 

S  m  1 

<j — 

6.33  x  107 

S  m1 

3  Parameters  are  obtained  from  Ref.  [7]  and  scaled  to  the  cell  size  in  this  article. 


i _ 

II 

1 

- 1 

_NZ_ 

/<rcos2a  +  /^sin2o: 

(kT  -  k^j  sinacosa 

0 


(k r  -  sinctcosce 

/<rsin2a  +  k^c  os2  a 

0 


rari 

ax 

ar 

ay 

_ 

SI3 

1 _ 

(33) 


therefore,  the  elements  of  K  can  be  obtained  from  Equation  (33): 


krcos2a  +  k^sin2a  (kT  -  k^j  sinacosa  0 
(kr  -  fc ) sinctcosct  /<rsin2a  +  kcC os2 a  0 


0 


0 


(34) 


2.9.  Solution  procedure 

A  flowchart  for  the  solution  procedure  using  our  coordinate 
transform  and  variable  extrusion  methods  is  presented  in  Fig.  11. 
The  charge  balance  equations  and  the  electrochemical  model  are 
solved  in  the  transformed  2-D  domain  for  the  electrode  pair,  and 
the  volumetric  heat  source  term  Qh  is  evaluated  as  a  function  of 
polar  angle  6  and  z  using  Equations  (19)— (25).  Then  Qjt\{6,z)  is 
transformed  as  a  function  of  Cartesian  coordinates,  Q^(x,y,z),  and 
extruded  to  the  3-D  jelly-roll  domain  using  the  algorithm  described 
by  Equations  (7)-(9).  The  thermal  model  is  coupled  with  the 
extruded  heat  source  term  and  solved  in  the  3-D  jelly-roll  domain. 
The  two  electrical  potentials,  &+  and  are  also  extruded  to  check 
the  accuracy  of  our  algorithms.  It  is  noted  that  the  variables 
extruded  from  the  2-D  domain  to  the  3-D  domain  must  be  explicit 
expressions  of  the  coordinates,  and  the  solutions  for  this  issue  are 
discussed  in  Appendix  B. 

The  geometries  for  the  two  computational  domains  are  pre¬ 
sented  in  Fig.  12,  the  3-D  domain  for  a  twenty-turn  jelly-roll  is 
regarded  as  a  continuous  phase  enclosed  by  two  spiral  surfaces 
following  the  polar  equation  (1).  For  the  interior  surface,  the  polar 
angle,  6,  changes  from  0  to  2tt;  for  the  exterior  surface,  6  changes 
from  (Nw  -  1  )2tc  to  Nw2tc,  where  Nw  =  20  is  the  number  of  turns  for 
the  wound  electrodes. 

Using  the  numerical  procedure  presented  in  Fig.  11,  only  Equa¬ 
tion  (26)  is  solved  in  a  3-D  domain  where  the  fluxes  are  continuous 
in  all  directions,  there  is  no  need  to  define  the  insulated  interior 
boundaries  in  the  3-D  domain,  and  the  meshing  requirements  are 
greatly  lowered.  Fig.  13  shows  the  mesh  pattern  for  the  3-D  jelly- 
roll  domain  as  our  method  is  applied,  only  21,950  domain  ele¬ 
ments  are  included,  and  the  computational  load  is  significantly 
reduced  compared  to  the  model  with  the  mesh  shown  in  Fig.  4. 


/ The  2-D  planar  surface 

a>  =o 

z\  \ 

1 — +$ 

i 

i 

Fig.  14.  Computational  domains  for  validating  the  coordinate  transform  algorithm. 


3.  Results  and  discussion 

The  mathematical  model  developed  in  Section  2  was  imple¬ 
mented  using  COMSOL  Multi-physics  4.3b  with  the  parameter 
values  listed  in  Table  2. 


3.1.  Validation  of  the  coordinate  transform  method 


To  validate  our  method  for  the  coordinate  transformation,  we 
need  to  show  that  Equation  (6)  solved  on  a  3-D  spiral  surface  and 
on  a  2-D  plane  have  the  same  solution.  In  this  case,  we  segregated 
the  charge  balance  equation  for  the  positive  electrode  current 
collector  by  setting  the  volumetric  current  source  as  a  constant 
which  is  equivalent  to  5C  discharge: 


(35) 


where  z'c  is  the  transversal  current  density  for  C  rate.  As  shown  in 
Fig.  14,  a  twenty-turn  spiral  surface  was  built  in  the  3-D  Cartesian 
coordinate  system;  a  rectangular  domain  was  built  in  the  ?-z  plane, 
the  height  of  the  rectangular  domain  is  same  with  that  of  the  spiral 
surface  and  the  span  of  the  rectangular  in  the  ?  direction  equals  the 
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▼  0 


Fig.  15.  Solutions  for  electrical  equation  solved  on  spiral  and  planar  surfaces. 


tangential  arc  length  of  the  spiral  surface.  The  boundary  conditions 
are  <P+  =  0  for  the  exterior  edge  of  the  spiral  surface  and  the  right 
edge  of  the  rectangular  domain.  With  a  constant  current  source,  the 
charge  balance  equation  is  steady-state.  The  results  from  the  spiral 
and  planar  surfaces  are  presented  in  Fig.  15,  and  there  is  perfect 
agreement  between  the  solutions  from  the  two  computational 
domains. 

3.2.  Validation  of  the  variable  extrusion  algorithm 

The  variable  extrusion  algorithm  (see  Equation  (10))  was 
validated  using  the  results  from  a  simulated  5C  discharge,  and 
the  purpose  was  to  check  agreement  between  a  variable  in  the  2- 
D  domain  and  its  extrusion  in  the  3-D  domain.  It  is  the  volu¬ 
metric  heat  source,  Qh,  that  couples  the  2-D  electrical  equations 
and  the  3-D  thermal  equations,  therefore,  the  Qh  values  evalu¬ 
ated  in  the  2-D  electrode  pair  domain  and  extruded  to  the  3-D 
jelly-roll  domain  at  the  beginning  and  end  of  5C  discharge  are 
presented  in  Fig.  16.  As  shown  in  Fig.  16,  the  heat  source  is  higher 
in  the  tab-contact  surfaces  of  the  two  electrodes  due  to  the  large 
current  density,  and  the  good  agreement  between  2-D  and  3-D 
values  confirms  the  high  accuracy  of  the  variable  extrusion 
algorithm. 

To  further  clarify  the  strength  of  this  algorithm,  the  Qh  vs  time 
plots  at  different  mapped  point  pairs  are  presented  in  Fig.  17.  A 
mapped  point  pair  includes  a  2-D  point  i  with  coordinates  (?,z) 
and  a  3-D  point  i"  with  coordinates  (x,y,z),  where  x,  y,  and  ? 
satisfy  the  relations  presented  in  Equations  (3)  and  (7)-(9).  The 


2-D  and  3-D  coordinates  for  the  three  mapped  point  pairs 
selected  at  different  coating  regions  are  listed  in  Table  3.  As 
shown  in  Fig.  17,  the  heat  source  values  in  the  NS-PD  (mapped 
point  pair  1)  and  ND-PS  (mapped  point  pair  3)  regions  are  higher 
than  those  in  the  ND-PD  (mapped  point  pair  2)  region;  the 
reason  is  that  the  NS-PD  and  ND-PS  regions  are  located  closer  to 
the  tab-contact  surfaces.  The  major  cause  for  the  error  between 
the  2-D  and  3-D  values  is  the  approximation  approach  presented 
in  Appendix  B,  and  the  accuracy  of  the  variable  extrusion  be¬ 
tween  the  2-D  and  3-D  domains  could  be  improved  by  including 
more  nodes  or  using  higher  order  basis  functions  in  the 
approximate  expressions  of  the  variables.  The  2-D  and  3-D  plots 
of  the  current  collector  potentials  for  the  positive  and  negative 
electrodes  at  the  end  of  5C  discharge  are  presented  in  Fig.  18;  and 
the  results  also  show  good  accuracy  in  the  2-D  to  3-D  variable 
extrusion. 

3.3.  Thermal  case  studies 

The  simulated  temperature  profiles  at  the  end  of  1C,  3C,  and 
5C  discharges  are  presented  in  Fig.  19,  where  the  initial  and 
ambient  temperatures  are  set  at  25  °C  and  the  convective  heat 
transfer  coefficient  at  the  exterior  spiral  surface,  ft,  is  set  as 
1.0  W  m-2  K-1.  As  shown  in  all  the  three  plots,  there  are  hot 
regions  located  at  the  outer  surface  and  near  the  end  of 
the  winding  turn,  which  are  corresponding  to  the  higher  heat 
generation  rates  in  the  tab-contact  surface  of  the  positive 
electrode. 
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(a)  Heat  generation  rate  at  the  beginning  of  5C  discharge 


(b)  Heat  generation  rate  at  the  end  of  5C  discharge 

Fig.  16.  Profiles  for  heat  generation  rates  (in  W  m-3)  in  2-D  and  3-D  domains  at  the  beginning  and  end  of  5C  discharge. 


The  effects  of  thermal  boundary  conditions  on  the  temperature 
profiles  are  also  discussed.  The  distributed  thermal  results  at  the 
end  of  5C  discharge  with  h  m  1.0  W  m-2  I<-1  and 
h  =  10.0  W  m-2  I<-1  are  presented  in  Fig.  20.  According  to  Fig.  20, 
the  temperature  of  the  jelly-roll  decreases  by  about  8  °C  when  the 


value  of  h  increases  from  1  to  10  W  m~2  K-1;  and  at  high  convection 
rate,  the  temperature  difference  between  the  center  and  outer 
surfaces  of  the  jelly-roll  becomes  larger  and  the  temperature  of  the 
outer  surface  drops  more  rapidly  than  the  center  surface. 

3.4.  Comparisons  between  our  method  and  conventional  method 

As  discussed  in  Section  2.2,  in  a  conventional  way,  the  electrical, 
electrochemical,  and  thermal  models  can  be  solved  in  an  extremely 
fine  meshed  3-D  domain  (Fig.  4).  However,  running  such  a  coupled 
multi-physics  model  exceeds  the  memory  limitations  for  most  PCs. 
For  this  reason,  we  present  comparisons  based  on  a  segregated 
thermal  model;  that  is,  the  heat  source  term  is  set  as  a  constant  (i.e. 
1  x  105  W  itT3)  and  thus  independent  of  the  electrochemical  and 


Table  3 

2-D  and  3-D  coordinates  for  different  mapped  point  pairs. 


x  [m] 

y  [m] 

z[m] 

f[m] 

Region 

Mapped  point  pair  1 

0 

0.025 

0.01 

0.01858 

NS-PD 

Mapped  point  pair  2 

0 

0.05 

0.03 

0.21307 

ND-PD 

Mapped  point  pair  3 

0 

0.082 

0.03 

0.63994 

ND-PS 
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A  2.8096 


▲  2.8096 


2.805 

2.8 

2.795 

2.79 

2.785 

2.78 

2.775 

2.77 

2.765 

2.764 


2.805 

2.8 

2.795 

2.79 

2.785 

2.78 

2.775 

2.77 

2.765 

2.764 


(a)  Electrical  potential  for  positive  electrode  current  collector 


(b)  Electrical  potential  for  negative  electrode  current  collector 

Fig.  18.  Profiles  for  electrical  potentials  (in  V)  of  current  collectors  in  2-D  and  3-D  domains  at  the  end  of  5C  discharge. 


1C  discharge 


I 


,  29.204 

29.202 

29.2 

29.198 
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29.192 

29.19 
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29.186 
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3C  discharge 


▲  39.853 

I  39.84 
■  39.82 

39.8 

39.78 

139.76 

39.74 
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5C  discharge 


▲  50.425 

I  50.4 
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50.3 
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L,s 

B50.I 

▼  50.099 


Fig.  19.  Profiles  for  jelly-roll  temperature  (in  °C)  at  the  end  of  1C,  3C,  and  5C  discharges. 
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electrical  models.  As  shown  in  Fig  21(a),  for  the  conventional 
method,  the  thermal  model  is  directly  solved  in  the  extremely-fine- 
mesh  domain  where  the  3-D  heat  source  term,  Q^(x,y,z),  is  set 
equal  to  constant.  As  shown  in  Fig.  21(b),  using  our  method,  the 
heat  source  term  is  defined  as  a  constant  in  the  2-D  electrode  pair 
domain  and  extruded  using  the  algorithm  developed  in  this  article 
to  a  coarsely  meshed  3-D  domain  where  the  thermal  model  is 
solved.  The  solutions  for  thermal  model  obtained  by  the  two 
methods  are  exactly  same  as  shown  in  Fig.  21.  By  the  conventional 
way,  it  takes  1  h  55  min  and  11.1GB  physical  memory  for  the  PC  to 
compute  the  results;  however,  with  our  method,  it  takes  only  3  min 
54  s  and  3.8  GB  physical  memory. 

If  the  electrical  and  electrochemical  models  are  included,  these 
additional  models  would  be  solved  in  the  extremely-fine-mesh  3-D 
domain  for  conventional  methods,  but  in  the  coarsely-meshed  2D 
Fig.  20.  Profiles  for  jelly-roll  temperature  (in  °c)  at  the  end  of  5C  discharges  with  domain  for  our  method.  Based  on  the  above  discussions,  we  can 
different  heat  transfer  coefficients. 


A  71.013 

I71 

■  70.98 
I  70.96 

70.94 

70.92 

70.9 

70.88 

■  70.86 
I  70.84 
I  70.82 

T  70.818 


(a)  Conventional  extremely-fine-mesh  method 


Temperature  at  t=1000s  (in  °C)  v  70.818 
(b)  Coordinate-transform  and  variable-extrusion  method 


Fig.  21.  Solving  a  segregated  thermal  model  using  (a)  conventional  extremely-fine-mesh  method  (b)  coordinate  transform  and  variable  extrusion  method. 
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conclude  that  our  method  would  be  advantageous  for  solving  the 
coupled  models. 

4.  Conclusion 

In  this  article,  a  numerical  method  is  presented  to  reduce  the 
computational  difficulties  in  modeling  a  spirally-wound  cylindri¬ 
cal  Li-ion  cell.  By  analyzing  the  geometric  relations  of  the  winding 
locus,  we  confirmed  that  the  charge  balance  equations  can  be 
solved  on  a  2-D  planar  surface  instead  of  a  3-D  spiral  surface 
without  losing  much  accuracy.  We  developed  a  coordinate  trans¬ 
form  and  variable  extrusion  algorithm;  and  with  these  the  elec¬ 
trical  and  thermal  variables  can  be  solved  and  evaluated  in  the  2-D 
electrode  pair  domain  and  then  delivered  to  the  3-D  jelly-roll 
domain  and  coupled  with  the  thermal  model.  The  accuracy  of 
variable  extrusion  could  be  improved  by  using  more  effective 
approximation  approaches  to  express  variables  as  explicit  func¬ 
tions  of  2-D  coordinates.  Our  numerical  method  enables  the  jelly- 
roll  model  to  be  solved  without  using  extremely  fine  mesh  and 
thus  greatly  lowers  the  memory  requirements  saves  computation 
time. 


Appendix  A 

Governing  equations  for  the  equivalent  circuit  model  (ECM) 

The  transient  response  equations  of  the  two  R—C  branches  in  the 
equivalent  circuit  presented  by  Fig.  7  are  as  follow: 


9Qp 

OQn 

at 


Qp 

ftpCp 


-El- 


On 

RnCn 


N,(0/) 


'  z'n,(0,0’) 


(Al) 


(A2) 


where  Qp  and  Qn  are  respectively  the  amount  of  charge  (in  the  unit 
of  C  m-2)  on  the  capacitors  Cp  and  Cn.  The  discharge  capacity  of  the 
cell  stack  (in  the  unit  of  Ah  m-2)  is  calculated  as  follow: 

QQdchg  _  lN,(fl/) 


dt  3600  [sh"1] 

The  depth-of-discharge  of  cell  (DOD)  is  defined  as: 

Qdchg 


DOD  = 


Qo 


(A3) 


(A4) 


where  Qo  is  the  theoretical  maximum  cell  capacity  per  unit  area  of 
electrode  pair.  The  potential  difference  over  the  cell  stack  is  related 
to  the  transversal  current  density  through  the  following  expression: 

&+(d)  -  <P-(d')  =  U  +  iN,(0/)^s  +  ^  (A5) 


Validation  of  the  ECM 

In  an  equivalent  circuit  model,  the  electrical  components  (Kp,  Cp, 
Rn,  Cn,  Rs,  and  U)  may  depend  on  the  depth  of  discharge  (DOD),  and 
these  DOD-dependency  relationships  can  be  empirically  deter¬ 
mined  by  fitting  the  ECM  to  real  data  or  results  from  a  full-order 
physics-based  model.  In  this  work,  the  ECM  was  validated  using 
the  pseudo-2D  porous  electrode  model  (P2D  model)  derived  by 
Doyle  et  al.  [15,16],  and  the  parameters  for  the  P2D  model  are  ob¬ 
tained  from  reference  7. 

The  open  circuit  voltage,  U ,  was  fit  to  a  C/100  discharge  profile 
generated  by  the  P2D  model  (as  shown  in  Fig.  A-l)  as  a  cubic 
spline  interpolate  function  of  DOD.  The  other  components  were  fit 
to  the  high  power  pulse  cycle  (EIPPC)  profiles  generated  by  the 
P2D  model  at  different  fixed  DOD  values.  The  HPPC  simulation 
protocols  include:  at  a  fixed  DOD,  a  5C  discharge  current  pulse  was 
first  applied  to  the  cell  for  10  s;  rest  the  cell  for  40  s;  a  3C  charge 
current  pulse  was  applied  for  10  s;  rest  the  cell  for  60  s;  and  then 
discharge  the  cell  at  1C  rate  for  to  next  DOD  value.  The  fits  of  ECM 
to  the  P2D-generated  EIPPC  profiles  at  a  specific  DOD  are  pre¬ 
sented  in  Fig.  A-2.  In  this  way,  the  estimated  values  for  electrical 
components  at  different  DOD  values  are  obtained  and  the  results 
are  listed  in  Table  A-l,  and  the  DOD  dependencies  of  these  com¬ 
ponents  on  DOD  are  also  approximated  using  the  cubic  spline 
interpolation. 


DOD 

Fig.  A-l.  Open-circuit  voltage  profile  generated  by  P2D  model. 


Time(s) 


Time(s) 


Fig.  A-2.  Fitting  ECM  to  the  HPPC  results  at  a  fixed  DOD  generated  by  P2D  model. 


Table  A-l 

Values  for  ECM  components  at  different  DOD 
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DOD 

Rs  [Q  m2] 

Rp  [fi  m2] 

Cp  [F  m-2] 

Rn  [Q  m2] 

C„  [F  m-2] 

1.524E-01 

2.150E-03 

1.058E-03 

3.307E+04 

2.804E-04 

6.242E+03 

2.580E-01 

2.150E-03 

1.127E-03 

2.837E+04 

2.828E-04 

5.383E+03 

3.635E-01 

2.152E-03 

9.919E-04 

3.277E+04 

2.272E-04 

6.678E+03 

4.691  E-01 

2.158E-03 

8.594E-04 

3.897E+04 

2.012E-04 

7.738E+03 

5.746E-01 

2.168E-03 

8.202E-04 

4.234E+04 

2.130E-04 

7.386E+03 

6.802E-01 

2.183E-03 

9.186E-04 

3.869E+04 

2.478E-04 

6.230E+03 

7.858E-01 

2.209E-03 

1.022E-03 

3.441  E+04 

2.821  E-04 

5.405E+03 

8.913E-01 

2.255E-03 

1.034E-03 

3.542E+04 

3.296E-04 

4.680E+03 

Appendix  B 

As  discussed  in  Section  2.9,  a  variable  extruded  from  the  2-D 
electrode  pair  domain  to  the  3-D  jelly-roll  domain  must  be  an 
explicit  expression  of  polar  angle  6  and  coordinate  z.  Let  W(?,z) 
stand  for  a  variable  continuously  distributed  in  the  2-D  domainf-z, 
and  *F(f,z)  can  be  approximated  as  the  following  expression: 


(bi) 

i  j 

where  W/j  is  the  value  of  W(J,z)  at  node  point  (?j,Zj),  and  p/(f )  and 
q/z)  are  linear  basis  functions  along  the  ?  and  z  directions.  The 
schematic  illustration  for  a  linear  basis  function,  p/(f ),  is  presented 
in  Fig.  B-l ;  according  to  the  plot,  pz(?)  =  1  at  node  point  ?  =  and 
Pi(?)  =  0  at  all  other  node  points  in  the  ?  direction;  and  in  the  two 
adjacent  intervals  of  node  point  ?  =  ?*,  and  [?i,?i+i],  p/(?) 

takes  linear  profiles.  Therefore,  p,-(?)  is  a  tent  function  that  can  be 
expressed  as  follow: 


(B2) 

According  to  Equation  (3),  £  is  an  explicit  function  of  6,  as  a 
result,  W(?,z)  can  thus  be  expressed  as  an  explicit  function  of  (9  and 
z. 


Fig.  B-l.  Plot  for  linear  basis  function. 
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List  of  symbols 

A+,  A_:  area  of  the  tab-contact  surfaces,  m2 

C:  rate  capacity  of  cell,  Ah 

Cp,  Cn:  capacitances  in  equivalent  circuit,  F  m-2 

Cp:  specific  heat  capacity  of  material,  J  K-1  kg-1 

D:  separation  distance  of  the  Archimedean  spiral,  m 

DOD:  depth  of  discharge 

h:  heat  transfer  coefficient  in  the  spaces  between  cells,  W  K-1  m-2 

io  transversal  current  density  for  C  rate,  A  m-2 

iN(0  0'):  transversal  current  density  through  a  cell  stack,  A  m-2 

iv>+,  volumetric  current  source,  A  m-3 

Iapp:  current  applied  to  or  from  the  module,  A 

k;  matrix  for  thermal  conductivities,  W  K-1  m-1 

kr:  thermal  conductivity  in  the  normal  direction  of  the  spiral  surface,  W  K-1  m-1 

kf  thermal  conductivity  in  the  tangential  direction  of  the  spiral  surface,  W  K-1  m-1 

nw:  index  of  spiral  winding  turns 

Nw:  number  of  winding  turns  in  a  jelly-rolll 

Nx,  Ny,  Nz:  heat  flux  in  the  x,  y,  and  z  directions,  W  m-2 

Nr:  heat  flux  in  the  normal  direction  of  the  spiral  surface,  W  m-2 

N%:  heat  flux  in  the  tangential  direction  of  the  spiral  surface,  W  m-2 

Qdchg-'  discharge  capacity,  Ah  m-2 

Qpc  e')  •’  volume-average  electrochemical  heat  source  through  a  cell  stack,  W  m-3 

Qjyffy  volume-average  Joule  heat  source  through  a  cell  stack,  W  m-3 

Qh-  Heat  generation  rate  in  the  electrode  pair  domain,  W  m-3 

Q^:  heat  generation  rate  extruded  to  the  jelly-roll  domain,  W  m-3 

Qp,  Qn :  amount  of  charge  on  capacitors,  C  m-2 

Rp,  Rn:  parallel  resistances  in  equivalent  circuit,  Q  m2 

r:  radius  coordinate,  m 

ro:  initial  radius  for  Archimedean  spiral,  m 

t:  time,  K 

T:  temperature,  s 

U:  open-circuit  voltage,  V 

x,  y,  z:  Cartesian  coordinates,  m 

a:  angle  between  the  normal  direction  of  spiral  surface  and  the  x  axis,  rad 

<5+,  5-:  thickness  of  current  collectors,  m 

5S:  thickness  of  cell  stack,  m 

<£>+,  O-:  electric  potentials  of  current  collectors,  V 

0:  polar  angle,  rad 

6e:  elevation  angle,  rad 

p:  density  of  cell,  kg  m-3 

<7_|_,  a  :  Electrical  conductivity  of  current  collectors,  S  m-1 
tangential  arc  length  of  the  spiral  surface,  m 


